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INTRODUCTION 


"Npmflpi 


In  an  earlier  report,  Ref.  [I],  the  present  authors  considered  the 
possibility  of  improving  the  tensile  modeling  of  geological  materials 
by  adopting  the  SRI  rate-dependent  brittle  fracture  model  BFRACT.  It 
was  concluded  that  the  model  was  too  costly  for  the  large-scale  ground 
shock  problems  of  interest,  due  to  the  large  number  of  memory  parameters 
required . 

It  was  recommended  then  that  versions  of  the  instantaneous  tensile 
fracture  model  commonly  in  use  be  re-evaluated.  These  versions  are 
mathematically  discontinuous  (since  some  stress  quantity  is  instantaneously 
reset  to  zero  when  the  tensile  strength  of  the  material  is  exceeded)  and 
hence  are  prone  to  stability-uniqueness  pioblems.  Such  procedures  have 
been  employed  for  many  years,  presumably  without  any  numerical  instabilities 
being  traced  to  their  use.  However,  it  has  been  noted  that  arbitrarily 
small  differences  in  input  data  do  lead,  on  occasion,  to  significantly 
altered  wave  forms.  It  was  further  recommended  that  a  plasticity  model 
be  tested  to  see  if  it  is  an  adequate  alternative  representation  for 
tensile  failure  in  ground  shock  calculations. 

To  study  these  recommendations  some  simple  one-dimensional  wave  pro¬ 
pagation  calculations  were  made.  The  anticipated  shortcomings,  for  both 
Instantaneous  fracture  and  plasticity  models,  are  documented.  No  method 
has  evolved  for  amending  the  rate-independent  models  to  overcome  these 
shortcomings . 
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Instead,  for  both  physical  and  theoretical  reasons,  a  rate-dependent 
model  for  brittle  tensile  fracture  in  ground  shock  calculations  is  proposed. 

The  model  adds  only  one  additional  state  variable,  the  relative  void 
volume  or  damage,  to  a  viscoplastic-like  description  of  material  in  tension. 

It  is  based  on  the  ductile  damage  accumulation  equations  advanced  at  SRI 
and  Sandia,  but  in  simplified  and  fully  incremental  form. 

Preliminary  uniaxial  strain  calculations  show  the  model  to  give  physically 
reasonable  results  and,  in  fact,  to  predict  stress-strain  paths  quite  close  to 
those  of  the  more  complicated  and  detailed  model  BFRACT.  Further 
testing  of  the  model,  for  more  realistic  two-dimensional  problems,  and  a 
study  of  cost  versus  benefit  in  ground  shock  calculations  is  in  progress. 
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II  RATE- INDEPENDENT  MODELS 

A.  Exact  Solutions 

A  compressive  step  pulse  of  width  k  and  amplitude  s  (negative  in 

compression)  is  assumed  to  propagate  toward  the  free  surface  of  a  semi-infinite 

bar  of  linear  elastic  material  (see  Fig.  1).  The  particle  velocity  v  is  then 

o 

(s/pc)  toward  the  left  with  p  the  undisturbed  density  and  c  the  wave  speed. 

The  exact  solution  of  the  problem,  for  times  greater  than  the  time  of  wave 
reflection,  is  critically  dependent  upon  the  strength  of  the  pulse  relative 
to  the  tensile  strength  of  the  material  Oq  . 

For  an  incoming  pulse  such  that  -s  <  O^a  tensile  wave  of  amplitude  -s 
reflects  back  down  the  bar. 

For  -s  >  C  the  ideal  solution  is  a  function  of  the  tensile  failure 
o 

option.  Assuming  instantaneous  fracture  with  all  stress  quantities  reset 
to  zero,  a  piece  of  length  &/2  spalls  off  with  uniform  velocity  equal  to 
twice  the  incoming  particle  velocity;  the  remainder  of  the  bar  is  at  rest. 

Assuming  plastic  tensile  failure  (as  implemented  in  the  cap  model  routine,  P.ef.  [2]), 
the  exact  solution  is  a  broadened,  reduced  tensile  pulse,  of  the  form  indicated  in 
Fig.  1,  returned  by  the  free  end  of  the  bar.  It  is  clear  that  the  plastic  solution 
is  not  a  physically  realistic  model  for  tensile  fracture.  (The  exact  solutions  for 
all  cases  under  consideration  are  illustrated  in  Fig.  1.) 

B .  Numerical  Tests 

Using  the  finite  difference  wave  propagation  code  WONDY,  Ref.  [3],  the 
numerical  solutions  obtained  for  the  step  pulse  problem  differ  from  the  exact 
solutions  discussed  above  in  ways  generally  characteristic  of  numerical 
solutions.  The  initial  square  wave  becomes  rounded.  Small  amplitude  numerical 
oscillations  appear,  causing  local  "overshoots"  and  "undershoots".  Where  a 

*)  Private  communication,  I.S.  Sandler,  Weidlinger  Associates. 
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Fig.  1  Step  Pulse  Propatnting  Along  a  Semi-Infinite  Bar  of  Linear  Elastic 
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Material 


spall  segment  forms,  its  velocity  profile  is  not  ideally  flat-topped. 

For  models  with  stability-uniqueness  problems,  these  numerical  deviations 
can  have  significant  consequences. 

Instantaneous  fracture  can  be  implemented  in  pure  continuum  fashion 

or  by  introducing  an  additional  node  and  permitting  separation  as  is  done  in  WONDY, 

when  the  tensile  strength  of  a  material  is  exceeded,  for  example.  Similar 

problems  arise  in  either  case.  The  former  option  has  been  selected  for  the 

present  numerical  test  along  with  the  following  non-dimensionalized  sample 

data:  p  =  1,  c  =  1,  0  =  0.0296,  £  =  0.4.  In  non-dimensionalized  units  then, 

the  initial  particle  velocity  is  equal  to  the  incoming  stress  pulse,  v  =  s. 

‘  o 

(Additional  details  of  the  code  implementation  are  given  in  Appendix  A.) 

*  ) 

In  Fig.  2  an  initial  square  pulse,  with  v  =  0.90  0  ,  is  shown. 

o  o 

At  a  later  time,  t  =  0.85,  the  particle  velocity  is  seen  to  have  approximately 
doubled  at  the  free  end  of  the  bar.  At  twice  the  previous  time,  t  =  1.7,  the 
wave  has  reflected  back  in  tension.  Ignoring  the  numerical  inaccuracies 
discussed  previously,  the  ideal  solution  for  a  pulse  of  amplitude  less  than 
the  tensile  strength  of  the  material  has  been  obtained. 

Increasing  the  initial  pulse  velocity  (synonomous  with  increasing  the 

stress)  a  small  amount  to  0.92  O  causes  the  ideal  solution  to  be  lost.  For 

o 

such  initial  data  the  third  time  frame  (t  =  1.7)  in  Fig.  3  shows  the  physically 
unrealistic  formation  of  a  spall  segment.  Instantaneous  tensile  fracture 
has  been  triggered  as  the  numerical  overshoots,  expected  to  develop  in  the 
finite  difference  solution,  produce  a  local  stress  level  in  excess  of  CT 


The  minus  sign  to  indicate  travel  toward  the  left  is  henceforth 
suppressed . 
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DISTANCE  ALONG  BAR 


3  Instantaneous  Tensile  Fracture.  v 


0.92  a 


A  small  change  in  input  data  is  seen  to  produce  a  large,  non-physical, 
change  in  the  results. 

Purely  computational  parameters,  unrelated  to  the  input  data  of  the 
physical  problem,  can  be  shown  to  unduly  influence  results  as  well.  These 
computational  parameters  arise  in  the  consideration  of  artificial  viscosity 
and  time  step  selection. 

Two  types  of  artificial  viscosity  are  added  in  WONDY  under  locally 
compressive  situations,  p  >  0.  The  quadratic  form  of  artificial  viscosity 

<1-  Bl<^*2 

is  most  effective  in  controlling  gradients  at  shocks  while  introducing 
mimimal  disturbances  elsewhere.  The  linear  form, 

q2  =  cp  Ax 

is  used  to  control  the  high  frequency  oscillations,  though  with  great  care 
as  this  form  is  prone  to  distort  a  general  solution  away  from  shock  areas 
as  well  as  close  in.  B^  and  B2  are  the  quadratic  and  linear  artificial  vis¬ 
cosity  coefficients  and  Ax  is  a  cell  width.  To  insure  general  numerical 
stability  the  time  step  used  to  advance  the  computation  is  controlled  by 
a  Courant-Friedrichs-Lewy  type  of  condition,  Ref.  [3], 


The  numerical  parameter  k  is  the  Courant  number,  lying  between  0  and  1. 

In  general  k  is  taken  close  to  1  to  maximize  the  possible  time  step,  while 
k  equal  to  1,  marginal  stability,  is  avoided.  Default  values  in  WONDY  for 
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B^,  B,,  and  k  are  2.0,  0.1  and  0.98,  respectively.  These  values  are  assumed 
unless  otherwise  stated. 

In  Fig.  4  the  incoming  pulse  velocity  is  again  v^  =  0.90  o  ,  but  k 
has  been  reduced  to  0.49  and  spall  once  more  appears.  Smaller  time  steps 
have  anomalously  produced  a  less  physically  accurate  result,  one  which  appears 
for  this  particular  configuration,  analogous  to  running  with  a  reduced 
amount  of  artificial  viscosity.  (This  is  not  a  general  result,  but  another 
manifestation  of  the  instantaneous  fracture  model's  theoretical  deficiency.) 

Spall  can  be  suppressed  (local  overshoots  reduced),  by  increasing  the 
strength  of  the  artificial  viscosity  while  all  other  quantities  are  held 
constant.  In  Fig.  5,  B2  has  been  increased  from  0.1  to  1.0,  vq  being 
0.90  0^  and  k  being  0.49.  With  this  set  of  parameters  a  tensile  wave  is 
once  again  observed  to  propagate  back  from  the  free  end  of  the  bar  as  required 
physically,  though  in  this  case  it  is  observably  damped. 

Table  1  summarizes  the  results  obtained  for  the  step  pulse  problem 
run  with  instantaneous  tensile  fracture  for  the  four  sets  of 
parameters  discussed  above. 


Table  1.  Step  Pulse  Problem  with  Instantaneous  Tensile  Fracture. 


Incoming 
Stress  Level 


Courant  Artificial  Viscosity  Spall  Fragment 
Number  Coefficient  Formation 

Quadratic  Linear 


DISTANCE  ALONG  BAR 


Fig.  4  Instantaneous  Tensile  Fracture.  vq  -  0.90  0 
Courant  Number  Halved. 


DISTANCE  ALONG  BAR 

Fig.  5  Instantaneous  Tensile  Fracture,  =  0.90  a Q 

Courant  Number  Halved.  Linear  Artificial  Viscosity 
Increased  Tenfold. 


When  the  square  wave  problem  is  solved  exactly  with  a  plasticity  model 

of  tensile  failure  no  spall  fragment  forms,  whatever  the  initial 

stress  level.  Since  the  set  of  constitutive  equations  is  well- 

behaved,  it  can  be  anticipated  that  the  numerical  results  will  closely  match 

the  theoretical.  At  an  incident  stress  level  50%  above  Oq  a  primary  tensile 

wave  of  width  i  and  amplitude  a  ,  followed  by  a  secondary  wave  of  width  £  and 

o 

/  *) 

amplitude  a  /2,  is  the  predicted  solution  .  This  wave  form,  as  shown  in 
Fig.  6,  is  approximately  obtained  numerically.  In  contrast,  the  required 
physical  solution,  spall,  is  seen  in  Fig.  7  where  instantaneous  tensile  fracture 
has  been  assumed  and  -s  is  again  50%  above  0  . 

As  expected,  for  this  highly  idealized  problem  (and  for  any  ground 
shock  calculation  where  spall  is  expected  to  be  important) , neither  rate- 
independent  model  of  tensile  failure  satisfies  both  the  requirements  of 
stability-uniqueness  and  of  physical  realism.  Rather  than  attempt  to 
amend  the  rate-independent  models,  the  choice  has  been  made  to  develop  a 
highly  simplified  rate-dependent  model  as  a  way  of  overcoming  the  discussed 
shortcomings. 


Refer  back  to  Fig.  1  for  the  general  form  of  the  plastic  failure 
solution  for  -s  >  a 

o 
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DISTANCE  ALONG  BAR 

Fig.  6  Plastic  Tensile  Failure,  v  =  1. 

o 


Instantaneous  Tensile  Fracture 


Ill  RATH- DEPENDENT  MODELS 


A.  Weidliiijjer  Visco-Damage  Model 

An  alternative  to  the  rate-independent  or  multi-parameter  rate- 
dependent  models  previously  discussed  is  suggested.  As  compared  to  SRI- 
BFRACT,  Ref.  [41,  or  to  the  brittle  spallation  theory  of  Davison  and 
Stevens,  Ref.  [5],  the  simplifying  assumption  is  made  that  a  single 
volumetric  state  variable  D,  the  tension-induced  relative  void  volume 
or  damage,  is  sufficient  to  characterize  tensile  fracture.  By  not 
recording  the  varying  sizes  and  orientations  of  planar  cracks  (actually 
found  in  bittle  materials)  it  is  felt  that  computational  costs  can  be 
sufficiently  reduced  to  make  a  rate-dependent,  and  hence  mathematically 
stable,  model  practical  for  ground  shock  calculat ions . 

It  is  proposed  that  the  incremental  elastic  constitutive  equation, 
o  =  be  modified  to  include  the  time-dependent  growth  of  damage  as 

foil ows : 

a  +  oj(a,D,c;t)D  =  M(o,D,r. ;  t)c ,  (1) 

where  a  is  the  average  stress  and  £  the  volumetric  strain.  The  damage 
growth  rate  D  is  taken  as  a  function  of  the  dependent  variables  and  time, 
but  not  of  their  rates 

D  =  f (o,D,t;t) 

To  satisfy  continuity  it  is  required  that  M  reduce  to  the  undamaged 
elastic  modulus  and  that  the  product  wD  go  to  zero  for  stresses  less 
tensile  than  some  threshold  value  CT and  for  damage  equal  to  zero.  The 
requirements  that,  in  a  uniaxial  strain  test  for  example,  time  to  fracture 
decrease  and  peak  tensile  strength  increase  witli  increasing  strain-rate 
impose  further  constraints  on  the  functional  forms  of  a,  D  and  M, 

(Such  rate  dependences  have  been  documented  for  rock  in  Ref.  [6],  for  example.) 


Unlike  the  various  nonlinear  u lastie-viseoplastie  constitutive  equations 


examined,  for  example  in  Refs.  [6,  7,  8],  this  form  has  been  tailored  to 
model  the  (rapid)  relaxation  of  stress  to  zero  characteristic  or  tensile 
failure,  provided  that  ;cD  remains  finite  (and  positive)  as  damage  increases. 

The  theory  of  spall  damage  accumulation  in  ductile  metals  of  Davison, 
Stevens  and  Kipp,  Ref.  [9],  of  Sandia  can  be  adapted  to  fit  the  above 
formalism.  The  Sandia  model  in  turn  relies  upon  the  SRI  fracture  models 
discussed  in  Ref.  [4],  for  spherical  void  size  distributions  and  growth 
laws  and  for  laboratory  data  to  set  the  material  parameters. 


The  time  rate  of  change  of  damage  separates  into  two  parts,  one 
dependent  on  void  nucleation,  the  other  on  void  growth 


D  =  D  +  D  =  f 
n  G 


(2) 


From  Ref.  [9| 


and 


D  =  c  v 
n  n*o 


0-0  +  o-o 

r  n  1  n  , 

exp[ - — - ]  -  1 


( 1-D) J 


(3) 


i>G  =  2  cG[2°  +  I0“°gI  "  la+0Gl]  D(1“D) 


where  is  the  initial  time  rate  of  change  of  void  volume,  the  average 

initial  (spherical)  void  volume,  0  the  nucleation  threshold,  O,  the 

n  1 

nucleation  sensitivity,  c  the  growth  coefficient  and  0  the  growth  threshold. 

Cr  G 

These  material  parameters  can  be  related  to  the  SRI  material  constants, 

T,  to  T£ ,  as  indicated  in  Appendix  B.  (Stress  and  strain  are  assumed 
1  o 

positive  in  tension.) 


An  integrated  form  has  been  assumed  in  Ref.  [9]  for  the  average  stress 
in  the  small  strain  approximation. 


a  =  Kc  (4) 

with  the  volumetric  strain  <:  and  the  bulk  modulus  K  both  taken  as  damage- 

dependent  quantities.  In  particular,  rCr  the  case  of  uniaxial  strain  and  for 

;ij,  a  reference  density, 

r(o,  D)  =  [(2  +  pr/p)(1-D)1/3  -  3] 
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has  been  obtained  .  The  bulk  modulus  of  a  composite  material  at  low  damage 

levels.  Ref.  [10],  is  assumed  with  V  the  Poisson  ratio  and  K  the  undamaged 

o  o 

modulus, 

3(l-v  ) 

K  ■  v  -  2<wVt  di 

o 

For  equal  to  1/3,  the  expression  becomes 

K  =  K  (1  -  3D)  (5) 

o 

In  this  case,  damage  is  constrained  to  be  less  than  or  equal  to  a  third. 

Differentiating  Eq.  (4)  with  respect  to  time  and  reverting  to  a  damage- 
independent  concept  of  strain,  one  may  write 


0  +  10  +  Kf 


ct  +  310  D  =  Kt;  ,  D  <  1/3  (6) 

The  following  identifications  can  then  be  made  for  the  unknown  functions 
in  Eq.  (1) 

u>(o,D,r.;t)  =  3K  e 
o 

M(o, D, f ; t )  =  K 

completing  the  specification  of  the  proposed  rate-dependent  volumetric 
stress-strain  relation  for  material  in  tension. 

Examining  Eqs.  (l)-(3),  it  is  clear  that  damage  initiation  is  controlled 

by  the  nucleation  threshold  a  for  the  usual  case  D(t=0)  =  0.  On  the  other 

n 

hand,  stress  relaxation  can  continue,  for  increasing  tensile  strain  and  large 
damage  values  (D  -»  1/3),  until  the  lower  of  the  two  threshold  stresses,  o 

n 

of  Op  is  reached.  Depending  upon  how  stringently  the  requirements  of  relaxation 

to  zero  stress  for  brittle  fracture  is  to  be  met,  a_  may  be  taken  anywhere  from 

G 

an  order  of  magnitude  or  two  less  than  a  to  zero,  to  give  the  appropriate 

n 


Actually,  the  Sandia  form  is  more  complex.  Thermal  effects  and  plastic 
flow  have  also  been  accounted  for  in  the  strain  components.  These  are 
being  disregarded  as  not  relevant  to  the  present  investigation. 
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asymptotic  stress  behavior  withouL  affecting  the  onset  of  damage 

*) 

accumulation . 

B .  SRI  Ductile  Fracture  Model 

While  the  routine  BFRACT  has  been  used  to  model  tensile  fracture  in 
a  ground  shock  test  calculation.  Ref.  [11],  the  closest  to  the  proposed 
rate-dependent  model  for  such  calculations  may  well  be  the  simpler, 
less  time  consuming  routine  DFRACT,  Ref.  [4],  constructed  for  ductile 
materials.  This  model  requires  the  calculation  of  solid  specific  volume 
and  pressure  in  addition  to  gross  specific  volume  and  pressure.  It 
accounts  for  the  progressive  growth  of  damage,  in  the  form  of  spherical 
voids,  and  for  the  relaxation  of  stress  under  the  influence  of  continuing 
tensile  strain.  Details  of  the  model  are  given  in  Appendix  B. 

The  main  advantage  of  the  proposed  Weidlinger  visco-damage  model  is 
that  it  is  expressed  in  incremental  form  and,  hence,  mathematically 
better  understood.  Also,  it  requires  the  storage  of  fewer  memory  parameters 
per  computational  cell,  corresponding  to  the  introduction  of  one, rather 
'•ban  two,  additional  state  variables. 

C .  Numerical  Tests 

If  a  constant  uniaxial  strain-rate  is  assumed,  e(t)  =  e  ,  the  set 

o 

of  Weidlinger  visco-damage  constitutive  equations  becomes 

dD  _  _f_ 

d£  ~  A 
°o 

¥  =  K  [1  -  3(D  +  #)] 
d£  o  f 

o 

with  the  function  f  given  by  Eqs.  (2)  and  (3)  and  D  in  the  range  0  to  1/3. 


The  parameter  cG  may  have  to  be  adjusted  as  well  to  achieve  full 
relaxation  within  a  reasonable  time  interval. 
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Adding  damage  accumulation  concepts  to  a  conventional  description  of 
the  deviator  stresses  limited  by  a  Mises  yield  condition,  as  in  Kef.  [4J, 
gives 


ds  .  =  2y  (dc .  .  -  i  6.  .  ^-V-) 
ij  ij  1  1.1  V 


2j2  =  £v,±!y2 


with 


y  =  U  (1  -  1.88D) 
o 

and 


Y  =  Y  (1  -  4D) 
o 

*) 

Numerically  integrating  the  resulting  set  of  equations  for  a  given 
material  yields  a  set  of  strain-rate  dependent  stress-strain  curves. 


The  results  of  such  numerical  integration  for  the  brittle  metal 

Armco  iron,  using  SRI  adapted  data  as  listed  in  Appendix  C,  are  shown  in 

Fig.  8.  These  axial  stress-strain  curves  conform  well  to  the  requirements 

set  forth  in  Section  A  for  brittle  tensile  fracture.  Stress  appears  to 

relax  fully  for  large  enough  tensile  strain  .  With  a  failure  strain 

inferred  from  Fig.  8  for  each  constant  strain  rate  curve,  time  to  fracture, 

Tj  =  t'.jA  ,  is  shown  in  Table  2  to  decrease  with  increasing  strain-rate, 

while  peak  tensile  stress  O  is  seen  to  increase  correspondingly. 

max 


*) 

Using  a  fourth  order  Runge-Kuita  scheme  with  the  Gill  modification, 
for  example. 

Jek  ) 

In  fact  stress  relaxes  to  0  a  low  value  of  stress  relative  to  ° 

G  max 

Alternatively,  to  get  true  relaxation  to  zero,  the  parameter  o„  can 
be  chosen  to  be  zero.  u 


18 


Table  2.  Uniaxial  Extension  Test,  Armeo  Iron,  Weidliuger  Visco-Damage  Model 


Strain  Rate 

o 

(l/sec) 

Time  to  Fracture 

Tf 

(.sec) 

Maximum  Stress 

q  9  max 

(10  dyn/cm“*0. 1  OPa ) 

io4 

,.7 

24. 

105 

0.  5 

54. 

io6 

<0.2 

71  . 

Though  limited  in  this  simple  numerical  test  to  uniaxial  strain,  the 
results  are  encouraging. 

The  stress-strain  curves  produced  by  BFRACT  driven  at  various 
constant  uniaxial  strain-rates  with  Armc.o  data  are  shown  superimposed  on  the 
visco-damage  curves  in  Fig.  9.  Corresponding  data  lias  been  used  Cor  both 

models  except  that  a  correction  has  been  made  to  A/  to  compensate  for  the 

ro 

change  from  plana r  to  spherical  voids.  (See  Appendix  B.)  The  earlier 
curves,  calculated  with  the  visc.o-damage  model,  are  in  excellent  agreement 
with  the  BFRACT  curves.  This  ,'.s  not  unexpected  since  in  a  one-dimensional 
strain  test  the  added  complexity  of  crack  directionality  goes  largely  unused 

As  a  second  trial  the  Weidlinger  visco-damage  model  and  BFRACT  have 
been  driven  at  various  strain-rates  to  simulate  the  uniaxial  extension  test 
with  data  representing  the  ductile  metal  aluminum.  Material  damage 
parameters  have  been  obtained  from  Ref.  [ 4 1  and  are  listed  in  Appendix  B. 

The  two  sets  of  curves,  shown  superimposed  in  Fig.  10,  are  in  qualitative 
agreement.  This  is  despite  the  fact  that  the  SRI  routine  DFRACT  would 
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TENS  STRESS 


SRI  FRACTURE  MODEL 
VI  SCO  -  DAMAGE  MODEL 


Fig,  g  Uniaxial  Strain  Rchnvior.  Armco  Iron 


*) 

have-  been  the  comparison  model  ot  choice  for  aluminum,  not  BFRACT 


Neither  set  of  curves  shows  the  sharp  drop  to  zero  stress  characteristic 

of  brittle  tensile  fracture,  due  to  the  choice  of  a  relatively  high  value  for 

a  ,,  (comparable  to  that  of  O  ).  Recall  from  Section  A  that  the  lower  of 
G  n 

the  two  threshold  values  controls  the  relaxation  of  stress  for  large  tensile 
strain  and  damage. 


*) 

The  distinction  between  planar  cracks  and  spherical  voids  has  been 

ignored  in  this  trial.  However,  for  aluminum,  the  models  are  not 

highly  sensitive  to  changes  in  the  value  of  the  fracture  parameter Y  . 

In  Fig.  11  Y  has  been  reduced  to  one-sixth  of  its  former  value, 
o 

all  other  quantitites  remaining  the  same.  Peak  tensile  stress  rises 
by  less  than  6%. 
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IV  CONCLUSION 


A.  Summary 

As  recommended  in  Ref.  [1],  the  rate-independent  models  -  instantaneous 
tensile  fracture  and  plastic  failure  -  have  been  re-examined  by  way  of  a 
simple  pulse  propagation  problem.  Neither  has  been  found  to  satisfy  both  the 
requirements  of  stability-uniqueness  and  of  physical  realism.  A  highly 
simplifed  rate-dependent  tensile  fracture  model,  utilizing  equations  and 
material  constants  from  SRI-DFRACT  and  the  Sandia  ductile  damage  accumulation 
model,  has  been  proposed  for  use  in  ground  shock  calculations.  It  is 
formulated  in  fully  incremental  form  and  requires  only  one  additional 
(volumetric)  state  variable.  Preliminary  uniaxial  strain  tests  have 
proven  encouraging. 

B.  Recommenda  t ions 

There  is  a  need  to  implement  and  test  the  Weidlinger  visco-damage 
model  in  more  realistic  two-dimensional  problems.  While  the  model  is 
expressed  in  terms  of  average  stress  and  volumetric  strain  and  is  thus 
generalizable  to  multi-dimensions,  it  is  important  to  test  the  hypothesis 
that  a  single  volumetric  variable  is  sufficient  to  significantly  improve 
the  modeling  of  brittle  tensile  failure.  That  it  improves  the  modeling, 
as  compared  to  the  rate-independent  models,  is  clear  since  it  is  mathemat¬ 
ically  well-behaved  and  yet  physically  more  realistic  than  plastic  tensile 
failure  (particularly  for  rocks).  Whether  the  additional  computer  cost  is 
justified  by  the  improvement  remains  to  be  further  documented. 

For  problems  in  which  crack  directionality  is  judged  to  be  especially 
important  it  may  be  possible  to  generalize  the  model  with  the  introduction 
of  an  additional  state  variable.  It  is  recommended  that  such  an  option  be 
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APPENDIX  A 


CODE  IMPLEMENTATION  OF  STEP  PULSE  PROBLEM 
Tension  Options 

The  step  pulse  problem  has  been  run  using  the  wave  code  WONDY  with  the 
usual  solid  equation  of  state  routine  STATl  simplified  to  model  a  linear  elastic 
material  and  the  three  tensile  failure  options  outlined  in  the  text: 
instantaneous  fracture  implemented  in  pure  continuum  fashion  and 
with  local  separation  of  the  mesh, and  plastic  tensile  failure.  There  are  three 
tension  model  options  available  in  WONDY.  The  first  (FCRIT  =  1) 
suppresses  fracture  and  allows  a  material  to  stretch  at  some  energy  dependent 
maximum  tensile  stress.  This  option  is  used  directly  to  model  plastic  failure 
with  the  maximum  tensile  stress,  a  constant,  equal  to  the  spall  strength 
of  the  material  CTq.  Option  1  can  also  be  used,  in  modified  form,  to  model 
instantaneous  fracture  in  a  continuum  manner.  In  this  case,  all  stress  quantities 
in  a  cell  are  reset  to  zero  when  the  stress  exceeds  the  spall  strength.  The 
second  WONDY  fracture  option  (FCRIT  =  2)  allows  two  zones  to  separate  whenever 

the  interpolated  stress  at  their  common  boundary  is  more  tensile  than  o^. 

,  The  third  WONDY 

This  is  the  non-continuum  approach  to  instantaneous  fracture, 
fracture  option  (FCRIT  =  3)  is  not  tested  in  the  present  study.  In  this  option, 

the  value  of  a  cumulative  damage  integral 

t 


K 


(a  -  a  )  dt, 


not  the  stress  itself,  is  compared  to  a  maximum  permitted  value  K 

max 

While  the  criterion  for  fracture,  K  >  K  ,  is  time  dependent,  the  implementation 

—  max 

of  fracture  in  terms  of  cell  boundary  separation  continues  to  be  instantaneous. 
Stress  is  not  relaxed  gradually  but  rather  instantaneously  reset  to  zero  when 
fracture  is  triggered.  Hence  option  3  can  be  expected  to  share  the  stability- 
uniqueness  problems  of  the  fully  time-independent  fracture  models  being  examined. 
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Non-d imensiona 1 Izat ion  and  Initial ization 

Since  there  are  no  dimensional  constants  in  WONDY,  running  in  non-dimensional 

form  is  suggested  for  parametric  studies.  A  < haracteristic  length  L  along 

the  bar  of  unit  value  is  chosen  as  the  reference  independent  variable.  The 

undisturbed  density  p  and  wave  speed  c  of  the  material  are  the  two  reference 

o  o 

dependent  variables  selected.  All  other  dependent  variables  and  time  are 
converted  to  their  appropriate  non-dimensional  forms. 


o  ->  o/p  c 
o  o 


2 


e  v- 

t  -> 


e/c  2 

o 

tc  /L 
o 


o 


where  c  is  the  energy  per  unit  mass  and  for  concreteness,  p  has  been  taken 

o 


as  2.7  gm/cm  ,  c  as  5.0  cm/sec,  and  a  as  2GPa. 

o  o 

The  mesh  is  set  up  with  100  uniform  cells  of  width  0.02.  The  step 
pulse  (of  width  0.40,  compressive  stress  amplitude  s,  and  density  1-s)  is 
initially  at  a  distance  0.6  from  the  free  end  of  the  bar  and  traveling  toward 
it,  as  illustrated  in  Fig.  A1  below. 

Free  Surfac 


o 

II 

> 

V  =  0 

Q 

II 

o 

Q 

II 

O 

p  =  1 

P  =  1— s 

p  =  l 

j=l  31  51  101 

Figure  Al.  Initialization  of  Step  Pulse  Problem 
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APPENDIX  B 


DFRACT  PRESSURE  VOLUME  RELATIONS 

The  equation-of-state  routine  DFRACT  is  based  on  a  damage  dependent 
set  of  pressure-volume  relations.  Damage  in  this  model  is  assumed  to 
take  the  form  of  spherical  voids  developing  under  the  influence  of  tensile 
pressure  (average  stress).  As  detailed  in  Ref.  [4], 


V  =  V  +  V 
s  v 


PV  =  P  V 
s  s 


P 

s 


ViT  -  +  T 

s  s 


dE  =  -PdV 


where  V  is  specific  gross  volume,  V  initial  specific  gross  volume,  V 

o  s 

specific  solid  volume  and  specific  void  volume.  The  pressure  associated 
with  the  gross  material  is  P,  that  with  the  solid  material  is  Pg.  The 
quantities  Kq,  F  and  E  have  their  usual  meanings  i.e.,  ambient  bulk  modulus, 
Gruneisen  parameter  and  internal  energy. 


The  time  rate  of  change  of  specific  void  volume  is  assumed  to  be 
caused  by  void  nucleation  and  void  growth 


V  =  V  +  V„ 
v  n  G 


where,  for  spherical  voids. 


V  H 
n 


3  P-T 

8ri3  T4  exP  “tT 


lo 


VG  ^ 


T. (P  -  T_)V 
I  2  v 


(P  has  been  assumed  positive  in  tension.) 


P  •>  t5 

p  <  t5 
P  >  T2 
p  <  t2 


The  relat ionships  between  the  SRI  material  parameters,  T  to  T,  , 

1  o 


and  the  constants  used  in  the  Weidlinger  visco-damage  model  are  listed 
in  Appendix  C.  Note  that  the  Grlineisen  parameter  can  be  set  to  zero  in  a 
given  trial  to  improve  the  correspondence  between  the  SRI  and  Weidlinger 


models. 


APPENDIX  C 


VISCO-DAMAGE  MATERIAL  PARAMETERS  FOR  ARMCO  IRON  AND  ALUMINUM 


Symbol  Relation  to  SRI  Units  Armco  Aluminum 

Parameter  Iron 


c 

n 

T4 

3 

No/cm  /sec 

12 

4.6x10 

9 

3x10* 

r° 

r  3 

8nT^ (spherical  voids) 

3 

era 

SxlO-13*) 

2.3x10' 

2 

_32(^-)T5T^  (planar 
cracks) 

E  =  Young's  modulus 

o 

n 

V  =  Poisson's  ratio 

T5 

dyn/cm2 

3xl09 

9 

3x10 

°i 

T6 

2 

dyn/cm 

4 . 6xl09 

4x10** 

CG 

Tj/3 

cm2/dyn/sec 

2xl0-4 

3.3x10' 

°G 

T2 

dyn/cm2 

2 . 0x1 08 

9 

4x10* 

*) 

SRI  data.  Ref.  [4],  has  been  used  directly  for  the  visco-damage  material 
parameters  except  in  the  case  of  y°  (Armco).  This  parameter  has  been 

modified  in  an  attempt  to  compensate  for  the  fact  that  the  nucleation 
size  parameter  T^  has  been  determined  for  planar  cracks,  but  is  being 

applied  to  spherical  voids  in  the  Weidlinger  model. 
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ATTN:  HNDED-SR 

U.S.  Army  Engineer  Div,  Ohio  River 
ATTN:  ORDAS-L 

U.S.  Army  Engr  Waterways  Exper  Station 
ATTN:  J.  Strange 
ATTN:  WESSA,  W.  Flathau 
ATTN:  J.  Zelasko 
ATTN:  WESSE,  L.  Ingram 
ATTN:  Library 
ATTN:  J.  Drake 
2  cy  ATTN:  WESSD,  J.  Jackson 

U.S.  Army  Material  S  Mechanics  Rsch  Ctr 
ATTN:  Technical  Library 

U.S.  Army  Materiel  Dev  &  Readiness  Cmd 
ATTN:  DRXAM-TL 

U.S.  Army  Missile  Command 
ATTN:  RSIC 

U.S.  Army  Nuclear  &  Chemical  Agency 
ATTN:  Library 

DEPARTMENT  OF  THE  NAVY 

Naval  Construction  Battalion  Center 

Civil  Engineering  Laboratory 
ATTN:  Code  L08A 
ATTN:  Code  L51,  R.  O'Dello 
ATTN:  Code  L51,  S.  Takahashi 

Naval  Electronic  Systems  Command 
ATTN:  PME  117-21 

Naval  Facilities  Engineering  Command 
ATTN:  Code  03T 
ATTN:  Code  04B 
ATTN:  Code  09M22C 

Naval  Material  Command 
ATTN:  MAT  08T-22 

Naval  Postgraduate  School 

ATTN:  Code  0142  Library 
ATTN:  G.  Lindsay 

Naval  Research  Laboratory 
ATTN:  Code  2627 

Naval  Surface  Weapons  Center 
ATTN:  Code  F-31 
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DEPARTMENT  OF  THE  NAVY  ^Continued) 

Naval  Surface  Weapons  Center 

ATTN:  Tech  Library  &  Info  Services  Br 

Naval  War  College 

ATTN:  Code  E-ll 

Naval  Weapons  Evaluation  Facility 
ATTN:  Code  10 

Officer  in  Charge 

Newport  Laboratory 

Naval  Underwater  Systems  Center 

ATTN:  Code  EM,  J.  Kalinowski 

Office  of  Naval  Research 
ATTN:  Code  715 
ATTN:  Code  474,  N.  Perrone 

Office  of  the  Chief  of  Naval  Operations 
ATTN:  OP  981 
ATTN:  OP  03EG 

Strategic  Systems  Project  Office 
Department  of  the  Navy 
ATTN:  NSP-43 

DEPARTMENT  OF  THE  AIR  FORCE 

Air  Force  Geophysics  Laboratory 
ATTN:  LWW,  K.  Thompson 

Air  Force  Institute  of  Technology 
ATTN:  Library 

Air  Force  Office  of  Scientific  Research 
ATTN:  J.  Allen 
ATTN:  W.  Best 

Headquarters 

Air  Force  Systems  Comnand 
ATTN:  DLW 

Air  Force  Weapons  Laboratory 
Air  Force  Systems  Cornua nd 

ATTN:  NTE,  M.  Plamondon 
ATTN:  NTES,  R.  Jolley 
ATTN:  NTES,  J.  Thomas 
ATTN:  NTES-C,  R.  Henny 
ATTN:  SUL 

Asst  Chief  of  Staff 
Intel  1 igence 

Department  of  the  Air  Force 
ATTN:  [NT 

Ballistic  Missile  Office 
Air  Force  Systems  Command 
ATTN:  MNNXH,  D.  Gage 
ATTN:  MNN 

ATTN:  MNNXH,  M.  Delvecchio 

Deputy  Chief  of  Staff 
Logistics  &  Engineering 
Department  of  the  Air  Force 
ATTN:  LEEE 


DEPARTMENT  OF  THE  AIR  FORCE  (Continued) 

Research,  Development  IS  Acq 
Department  of  the  Air  Force 
ATTN:  AFRDQSM 

Foreign  Technology  Division 
Air  Force  Systems  Comnand 
ATTN:  tills  Library 

Rome  Air  Development  Center 
Air  Force  Systems  Comnand 
ATTN:  TSLD 

Strategic  Air  Command/XPFS 
Department  of  the  Air  Force 

ATTN:  NRI-STINFO  Library 

VELA  Seismology  Center 
ATTN:  G.  Ullrich 

DEPARTMENT  OF  ENERGY 

Department  of  Energy 
Albuquerque  Operations  Office 
ATTN:  CTID 

Department  of  Energy 
Nevada  Operations  Office 

ATTN:  Mail  &  Records  for  Tech  Library 

DEPARTMENT  OF  ENERGY  CONTRACTORS 

Lawrence  Livermore  National  Laboratory 
ATTN:  Tech  Info  Dept  Library 

Los  Alamos  National  Scientific  Lab 
ATTN:  MS  364 
ATTN:  R.  Bri dwell 
ATTN:  MS  670,  J.  Hopkins 

Oak  Ridge  National  Laboratory 
ATTN:  Civil  Def  Res  Proj 
ATTN:  Central  Research  Library 

Sandia  National  Labs 
ATTN:  3141 
ATTN:  A.  Chabai 
ATTN:  L.  Hill 

Sandia  National  Labs 
Livermore  Lab 

ATTN:  Library  &  Security  Class  Div 

OTHER  GOVERNMENT  AGENCIES 

Central  Intel! igence  Agency 
ATTN:  OSWR/NED 

Department  of  the  Interior 
Bureau  of  Mines 

ATTN:  Tech  Library 

Federal  Emergency  Management  Agency 

ATTN:  Hazard  Eval  &  Vul  Red  Div 
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DEPARTML'i'  OF  DEFENSE  CONTRACTORS 


DEPARTMENT  OF  DEFENSE  CONTRACTORS  (Continued) 


Aerospace  Corp 

ATTN:  Technical  Info  Services 


Institute  for  Defense  Analyses 
ATTN:  Classified  Library 


Agbabian  Associates 

ATTN:  M.  Agbabian 


J.  H.  Wiggins  Company,  Inc 
ATTN:  J.  Collins 


Applied  Theory,  Inc 
2  cy  ATTN  J.  Trul io 


AVCO  Research  8  Systems  Group 

ATTN 

Library  A830 

BDM  Corp 

ATTN 

T.  Neighbors 

ATTN 

Corporate  Library 

Boeing  Company 

ATTN 

Aerospace  Library 

Cal ifornia 

Institute  of  Tech 

ATTN 

T.  Ahrens 

Cal ifornia 

Research  8  Technology 

ATTN 

M.  Rosenblatt 

ATTN 

S,  Schuster 

ATTN 

Library 

ATTN 

K.  Kreyenhagen 

California  Research  8  Tech,  Inc 
ATTN:  0.  Orphal 


Cal  span  Corp 

ATTN:  Library 

University  of  Denver 

Colorado  Seminary 

ATTN:  J.  Wisotski 

EG8G  Washington  Analytical  Srv  Center,  Inc 
ATTN:  Library 

Eric  H.  Wang 

Civil  Engineering  Rsch  Fac 
ATTN:  N.  Baum 

Ga rd ,  Inc 

ATTN:  G.  Neidhardt 

General  Electric  Company— TEMPO 
ATTN:  DAS  I  AC 

Higgins,  Auld  8  Associates 
ATTN:  N.  Higgins 
ATTN:  H.  Auld 
ATTN:  J.  Bratton 

Higgins,  Auld  8  Associates,  Inc 
ATTN:  S.  Blouin 

1 1 T  Research  Institute 

ATTN:  Documents  Library 
ATTN:  M.  Johnson 
ATTN:  R.  Welch 


Kaman  AviDyne 

ATTN:  N.  Hobbs 
ATTN:  Library 

Kaman  Sciences  Corp 
ATTN:  Library 

Lockheed  Missiles  8  Space  Co,  Inc 
ATTN:  T.  Geers 
ATTN:  Tech  Info  Center 

Lovelace  Biomedical  8  Env  Rsch  Inst,  Inc 
ATTN:  R.  Jones 

McDonald  Douglas  Corp 

ATTN:  R.  Halprin 

Merritt  CASES,  Inc 

ATTN:  J.  Merritt 
ATTN:  Library 

Nathan  M.  Newmark  Consult  Eng  Svcs 
ATTN:  N.  Newmark 

Physics  International  Company 
ATTN:  Technical  Library 
ATTN :  E .  Moore 
ATTN:  L.  Behrmann 
ATTN:  F.  Sauer 
ATTN:  J.  Thomsen 

R  &  D  Associates 

ATTN:  J.  Lewis 
ATTN:  J.  Carpenter 
ATTN:  W.  Wright,  Jr 
ATTN:  R.  Port 
ATTN:  C.  MacDonald 
ATTN:  Tech  Info  Center 
ATTN:  P.  Haas 

Science  Appl ications ,  Inc 

ATTN:  Technical  Library 

Science  Applications,  Inc 
ATTN:  D.  Maxwell 
ATTN:  D.  Bernstein 

Science  Appl ications ,  Inc 
ATTN:  W.  Lay son 

Southwest  Research  Institute 
ATTN:  W.  Baker 
ATTN:  A.  Wenzel 

Systems,  Science  &  Software,  Inc 
ATTN:  J.  Murphy 
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DEPARTMENT  OF  DEFENSE  CONTRACTORS  (Continued) 
SRI  International 


ATTN 

ATTN 

ATTN 

ATTN 


Y.  Gupta 
G.  Abrahaiuson 
ii.  Gasten 
D.  Reouyh 


s> steins,  Science  &  Software,  Inc 


ATTN 

T. 

R  i  ney 

ATTN 

D. 

Grine 

ATTN 

L ibrary 

ATTN 

T. 

Cherry 

Tek, 

Inc 

ATTN 

. 

Green 

ATTN 

A. 

Abou-S< 

ATTN 

Lit 

;rary 

Teen 

Inc 

ATTN 

L. 

Hwang 

ATTN 

Library 

DEPARTMENT  OF  DEFENSE  CONTRACTORS  (Continued) 
TRW  Defense  &  Space  Sys  Group 


ATTN 
ATTN 
2  cy  ATTN 


P.  Bhutta 
Tech  Information  Center 
N.  Lipner 


TRW  Defense  &  Space  Sys  Group 
ATTN:  P.  Dai 
ATTN :  £ .  Wong 

Universal  Analytics,  Inc 
ATTN:  E.  Field 

Weidlinger  Assoc,  Consulting  Engineers 
ATTN:  M.  Baron 

Weidlinger  Assoc,  Consulting  Engineers 
ATTN:  J.  Isenberg 


[ 


